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Broken spatial inversion symmetry in spin-orbital coupled systems leads to a mixing between orbitals with different 
parity, which results in unusual electronic structures and transport properties. We theoretically investigate the possibility 
of multipole ordering induced by a parity mixing. In particular, we focus on the system in which the parity mixing 
appears in a sublattice-dependent form. Starting from the periodic Anderson model with such a local parity mixing, we 
derive an extended Kondo lattice model with sublattice-dependent antisymmetric exchange couplings between itinerant 
electrons and localized spins. By the variational calculation, simulated annealing, and Monte Carlo simulation, we show 
that the model on a quasi-one-dimensional zig-zag lattice exhibits an odd-parity multipole order composed of magnetic 
toroidal and quadrupole components at and near half filling. The multipole order causes a band deformation with the 
band bottom shift and a magnetoelectric response. The results suggest that unusual odd-parity multipole orders will be 
widely observed in multi-orbital systems with local parity mixing. 


1. Introduction 

The antisymmetric spin-orbit coupling, which is induced 
by inversion symmetry breaking of the crystal structure, has 
attracted interest because it leads to fascinating phenomena, 
such as the magnetoelectric effect'^^ and noncentrosymmetric 
superconductivity.^ " In such noncentrosymmetric systems, 
the Hamiltonian for the antisymmetric spin-orbit coupling is 
written by 

‘Hasoc(*) = ■ o-~ X VVpot) ■ rr, (1) 

where gik) is the polar vector with respect to the wave vector 
k, cr is the spin operator, and Wpot is a potential gradient 
due to the inversion symmetry breaking. This antisymmetric 
spin-orbit coupling often provides a clue for understanding of 
peculiar electronic and transport properties, such as the spin 
splitting in the band structure. 

In such spin-orbital coupled phenomena, not only uniform 
but also local (site) inversion symmetry breaking has recently 
drawn considerable attention.*^The local inversion sym¬ 
metry breaking means that the inversion symmetry is broken 
at every lattice site in a different manner from site to site. For 
instance, a two-dimensional honeycomb lattice breaks the lo¬ 
cal inversion symmetry at every lattice site in an alternating 
way on the two sublattices, while it preserves the inversion 
symmetry at each center of the nearest-neighbor bonds and 
the hexagons. In these systems, the parity mixing between dif¬ 
ferent orbitals occurs locally in a sublattice-dependent form, 
which we call the local parity mixing. In the presence of the 
local parity mixing, g{k) in Eq. ( 1 ), also depends on the sub¬ 
lattice. This suggests that new spin-orbital coupled phenom¬ 
ena are expected by using the sublattice degree of freedom 
in crystals. For example, once a sublattice-dependent mag¬ 
netic or electronic order is induced in the system with local 


parity mixing, multipoles extending over different sublattice 
sites are simultaneously activated in odd-parity sectors, such 
as a magnetic quadrupole and an electric octupole.'*’ 

Such unusual odd-parity multipoles have been recently 
studied for unconventional superconductivity^°“^^ and elec¬ 
tronic orders.A typical case was theoretically dis¬ 
cussed on a zig-zag chain, where the local inversion sym¬ 
metry is broken and the sublattice-dependent antisymmetric 
spin-orbit coupling is present.*^ In this case, an antiferromag¬ 
netic order can be regarded as an odd-parity multipole order, 
accompanying both magnetic toroidal and quadrupole com¬ 
ponents.'*'® Accordingly, the antiferromagnetic order ac¬ 
quires an unusual magnetoelectric coupling to an electric cur¬ 
rent.Therefore, such odd-parity multipoles are intriguing 
for exploring new spin-orbital coupled phenomena. However, 
the microscopic theory is not fully investigated. In particular, 
it is still unclear when and how such multipole orders are sta¬ 
bilized. It is highly desired to develop a general framework 
of the microscopic theory for the systems with local parity 
mixing. 

Experimentally, the local inversion symmetry breaking is 
widely found in multi-orbital systems on particular lattice 
structures. For instance, there are several /-electron materi¬ 
als possessing the lattice structures without local inversion 
symmetry, such as ferromagnetic superconductors UGe2,^* 
URhGe,*'** and UCoGe,*''^*® the zig-zag chain compounds 
LnM2A\\Q (Ln-Ce, Nd, Gd, Dy, Ho, and Er; M=Fe, Ru, 
and Os),*®^ the distorted honeycomb compounds a- and /- 
YbAlB4,‘'^^^ and the diamond-structure compounds RT2X20 
(R=Pr, La, Yb, and U, T^Fe, Co, Ti, V, Nb, Ru, Rh, and 
Ir, and X =A 1 and Zn).^*“^® Although these materials can be 
candidates for unusual odd-parity multipole ordering through 
the local parity mixing, their properties have not been studied 
from such a viewpoint. In order to stimulate experiments and 
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theories, it is important to clarify the microscopic mechanism 
for odd-parity multipole ordering. 

In this paper, we investigate a microscopic model by taking 
into account the hybridization between conduction and local¬ 
ized orbitals with different parity; e.g., between a d compo¬ 
nent in conduction electrons and localized / orbitals. Starting 
from the periodic Anderson model with the antisymmetric hy¬ 
bridization defined on a quasi-one-dimensional lattice struc¬ 
ture composed of zig-zag chains, we derive an effective low- 
energy model with antisymmetric exchange couplings be¬ 
tween conduction electrons and localized spins. This is an ex¬ 
tension of the Kondo lattice model to the case with local parity 
mixing. We show that the effective antisymmetric exchange 
couplings stabilize a Neel-type antiferromagnetic order with 
the spin polarization perpendicular to the zig-zag plane. This 
is an odd-parity multipole order, which gives rise to a pecu¬ 
liar band deformation and magnetoelectric effects.Using 
several numerical calculations, we clarify that the multipole 
order is widely stabilized at and near half filling in the ex¬ 
tended Kondo lattice model. 

The organization of this paper is as follows. In Sect. 2, we 
describe the derivation of the Kondo lattice model including 
antisymmetric exchange couplings from an extended periodic 
Anderson model. In Sect. 3, we clarify how the antisymmetric 
exchange couplings favor a multipole order associated with 
an antiferromagnetic order. We also compute the electronic 
band structure and magnetoelectric effect under this multi¬ 
pole ordering. Moreover, we examine the stability of the mul¬ 
tipole order systematically by the variational calculation for 
the ground state, simulated annealing, and Monte Carlo sim¬ 
ulation at finite temperatures. The last section is devoted to a 
summary of this paper. The canonical transformation leading 
to the extended Kondo model is given in Appendix. 

2. Model 

2. 1 Periodic Anderson model with antisymmetric hybridiza¬ 
tion 

In this paper, we consider the effect of the antisymmetric 
spin-orbit coupling in the presence of local parity mixing. 
For this purpose, we introduce an extended periodic Anderson 
model with the antisymmetric hybridization between different 
parity orbitals,whose Hamiltonian is given by 

'K = 'Ko+'Hi, (2) 

where 

■Ho = - ^{tijc]^Cj^ + H.c.) + U ^ ^ (3) 

i i,(T 

■Hi = Vi(k)(clj„^ + H.c.) + gf(k) ■ sf(k). (4) 

l,k,(T l,k,(T 

Here, (C;^) and /.^ (f.^) are the creation (annihilation) op¬ 
erators of conduction and localized electrons with spin cr at 
site i = (p, 1) {p and I denote the indices for the unit cell 
and sublattice, respectively); and are their 

Fourier transformations, respectively. The first term in Eq. (3) 


represents the kinetic energy of conduction electrons; tij is 
the hopping matrix element from site j to i. The second and 
third terms are the on-site Coulomb interaction for localized 
electrons and the atomic energy of localized electrons, re¬ 
spectively Meanwhile, the first term in Eq. (4) 

represents the (symmetric) hybridization between conduction 
and localized electrons. We note that Ho and the first term in 
Hi comprise the conventional periodic Anderson model.®' 

On the other hand, the second term in Eq. (4) is introduced 
to describe an antisymmetric hybridization between c and / 
electrons; g^/{k) is the polar vector describing the antisym¬ 
metric part, and 

+ H.C.), (5) 

cr,(T' 

where cr = (cr^, cr^, cr^) is the vector of Pauli matrices, as in 
Eq. (1). This term originates from the cooperation between 
the odd-parity crystalline electric field, atomic spin-orbit cou¬ 
pling, and off-site hybridization between orbitals with differ¬ 
ent parity.Although similar hybridizations may appear 
between the same c-c or /-/ orbitals, we omit them and con¬ 
sider only the hybridization between different orbitals c and 
/ in the present model. This is justified in some realistic situ¬ 
ations, e.g., when the atomic spin-orbit coupling for conduc¬ 
tion electrons is weak and the intersite overlap integrals be¬ 
tween localized electrons are small. 

The model in Eq. (2) is applicable to the generic systems 
without spatial inversion symmetry. In the present study, we 
focus on a specific situation where the antisymmetric hy¬ 
bridization originates from local parity mixing. As mentioned 
in the introduction, the local parity mixing exists in the lat¬ 
tice structures in which the inversion symmetry is broken 
at each site in a sublattice-dependent form; typical exam¬ 
ples are a zig-zag chain, honeycomb lattice, and diamond lat¬ 
tice. On these lattices, the sublattice-dependent antisymmet¬ 
ric hybridization is present because each site is affected by the 
sublattice-dependent odd-parity crystalline electric field. 

In the following sections, we consider one of the simplest 
realizations of the local parity mixing, a three-dimensional 
system composed of weakly-coupled one-dimensional zig¬ 
zag chains, as shown in Eig. 1. We set the x and y axes in 
the plane on which the chains lie, with taking the x axis in the 
chain direction; we set the lattice constants a = b = c - 1. 
Eor this setting, the inversion symmetry is broken at each lat¬ 
tice site, and the odd-parity crystalline electric field is present 
along the y direction with an alternating sign for two sublat¬ 
tices. Then, the local parity mixing results in the sublattice- 
dependent antisymmetric vector (k). Eor the current sys- 

C f 

tern with the quasi-one-dimensional structure, g^ (k) is ap¬ 
proximately given by 

gf(k)^(-l)‘gzsmkjc, (6) 

where I = 0(1) for the A(B) sublattice, g is the coupling con¬ 
stant, and z is the unit vector in the z direction. The form of 
gi (k) II Z is understood from Eq. (1) by considering the quasi 
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Fig. 1. Schematic pictures of (a) a single zig-zag chain and (b) a three- 
dimensional lattice composed of the zig-zag chains. The zig-zag chain is in 
the xy plane: the x and y axes are taken in the parallel and perpendicular direc¬ 
tions to the chain, respectively. The dashed rectangle represents the unit cell 
including two sublattice sites, A and B. a, b, and c are the lattice constants; 
t[, t2, ?3, and ^4 represent the hopping integrals. 


one-dimensional kinetic motion k\\ x and VVpot || y. Note that 

C f 

gj (k) acts only between c and / electrons on the same sub¬ 
lattice. 

2.2 Extended Kondo lattice model 

With these preliminaries in the previous subsection, we de¬ 
rive an effective low-energy model by treating Elx as the per¬ 
turbation to "Kq. This is done by a standard procedure on the 
basis of the canonical transformation, that is the Schrieffer- 
Wolff transformation.®^ Note that the procedure for the model 
without the antisymmetric hybridization term (i.e., g = 0) 
leads to the standard Kondo lattice model.®^^® In other words, 
we here extend the Kondo lattice model by including the ef¬ 
fect of the antisymmetric hybridization. In the derivation, we 
assume that the number of / electrons is one at each site, as in 
the standard procedure. We also assume that the hybridization 
Vi{k) has only the on-site component and sublattice indepen¬ 
dent; we drop the k and I dependences. By the straightforward 
calculations given in Appendix in detail, we obtain the effec¬ 
tive Hamiltonian as 

'Hjx-KLM = + H.C.) + — ^ cl^CTo-o-'Cjg., ■ Si 

ij,cr i^cT,cr' 

+ Yj {T^kySliS]^.^-S-iS^^,^ + Slinik>k)6R^+n.c\ 
l,p,k,k' ^ 

+ -Spisjk,k - Spi4k,k')^R^, 

hp,k,k' 

(7) 

where s, = sink, 6r^ = ^m'k = 

^^Ik'k ~ ~ 

~ represents the localized spins at 

the unit cell p and sublattice I and 5*^ = 5^^ + i5^,. 


In Eq. (7), the first and second terms are the kinetic energy 
of conduction electrons and the on-site exchange coupling be¬ 
tween localized spins and conduction electrons, respectively. 
These two terms represent the standard Kondo lattice model. 
Here, J is in the second order of V, as shown in Eqs. (A ll) 
and (A-14). In the first term, we take into account four dif¬ 
ferent hopping matrix elements; intrachain hoppings between 
nearest and next-nearest neighbor sites, fi and t 2 , respectively 
[see Fig. 1(a)], and interchain hoppings between the same 
sublattices for the neighboring sites in the y and z directions, 
f 3 and f 4 , respectively [see Fig. 1(b)]. 

The third and fourth terms in Eq. (7) are the antisymmet¬ 
ric exchange couplings, derived from the antisymmetric hy¬ 
bridization in Eq. (4). The third term results from the combi¬ 
nation of V and g^k). The coefficient Di is given by 

A = (-!)'VTg, (8) 

where G is proportional to [see Eqs. (ATS) and (AT6)]; 
Di is proportional to gV. The st^ - sink^ dependence in 
Eq. (7) comes from Eq. (6), which is the source of the band 
deformation and magnetoelectric effect, similar to the case 
of toroidal ordering;this is demonstrated in Sect. 3.2. 
The fourth term in Eq. (7) comes from the second order of 
g‘i^{k). It plays no essential role in the peculiar electronic 
and transport properties, since it is symmetric with respect to 
(k, k) —> (-k, -k'). Hereafter, assuming the situation V > g, 
we mainly consider the case with J > |D;| > G; we take J and 
G to be positive. For simplicity, we treat S, as a classical spin 
with |S,| = 1/2 in the following analysis. 


3. Results 

3.1 Two-site problem 

Let us hrst examine the effect of the third and fourth terms 
in Eq. (7) by considering a simple two-site problem. We dis¬ 
cuss what type of spin configuration is stabilized by these an¬ 
tisymmetric exchange couplings arising from the antisymmet¬ 
ric hybridization. We consider two sites connected by t 2 in a 
single chain, 1 and 2 (belonging to the same sublattice, say 
A-sublattice), as shown in Fig. 2(a). For the two sites, the ex¬ 
change couplings in Eq. (7) are explicitly written in the matrix 
form as 


where 


qxn _ 
'”2 site “ 


_ 

'”2 site “ 


7^2 site ~ 


_ 2 site 

P(12 \ 

■"2 site 

2 site 

^22 h 

■"2 site / 

cr^ H- 1 

,2 ' 2 


czs 0 .VTG 

s^45^ +S\)a- + 1 


(9) 

( 10 ) 

-s-Dcjy. (11) 


“^ILe = ”^24 and is given by -//^"ite with 1 ^ 2. 

We take the basis (ci-f, ci^, C 2 'f, C 2 x) in Eq. (9), and cr® is the 
2x2 unit matrix. Here, we omit the y component of localized 
spins, Sy, without loss of generality, as it gives an equivalent 
contribution to that from due to the spin rotational sym- 
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metry in the xy plane. Equation (9) indicates that there is an 
effective hopping of conduction electrons, induced by the an¬ 
tisymmetric exchange coupling proportional to V7g. The ef¬ 
fective hopping depends on the directions of localized spins, 
as shown in Eq. (11). This differentiates the energies for dif¬ 
ferent spin configurations. 



\/\/\/\ 


Fig. 2. Schematic pictures of magnetic orders on the zig-zag chain: (a) the 
antiferromagnetic order with up-down spin moments along the z direction 
(z-UD) and (b) the up-up-down-down antiferromagnetic order with the mo¬ 
ments along the x direction (jr-UUDD). In (a), the dashed oval shows the two 
neighboring sites in the same sublattice, 1 and 2 . 


We first assume two different types of ferromagnetic con¬ 
figurations for the two A-sublattice sites: one is the state with 
spin polarization in the x direction, and the other in the z direc¬ 
tion. Eor these two cases, the matrix elements of the Hamilto¬ 
nian in Eq. (9) are given by 



The eigenvalues are obtained as 





(14) 

(15) 


The results in Eqs. (14) and (15) clearly show that the an¬ 
tisymmetric exchange couplings prefer to the ferromagnetic 
configuration with magnetic moments in the z direction for 
the two A-sublattice sites. 

Next, let us consider a Neel-type antiferromagnetic spin 
configuration for the A-sublattice sites. In this case, the 
Hamiltonian is given by 




for the cases with the antiferromagnetic moments along the x 
and z directions, respectively. The eigenvalues are obtained as 

(^x) = + ( ^ -H ^ S^ J + 2^ Si,, (18) 

( 19 ) 

Thus, the stability of the x and z orders is opposite to the fer¬ 
romagnetic case: the antisymmetric exchange couplings favor 
the antiferromagnetic configuration with magnetic moments 
in the x direction. 

On the other hand, the spin configuration between the A 
and B sublattices will depend on an effective exchange cou¬ 
pling induced by the intersublattice hopping fj. Eor instance, 
when fi is much larger than fa, the antiferromagnetic config¬ 
uration between two sublattices is induced at and near half 
filling in the strongly-correlated region. This is partly under¬ 
stood by considering the second-order perturbation with re¬ 
spect to fi at half filling, which leads to an effective antiferro¬ 
magnetic interaction between localized spins. Meanwhile, the 
ferromagnetic configuration is favored at low and high fillings 
due to the double-exchange mechanism, which is an effective 
ferromagnetic interaction between localized spins induced by 
the kinetic motion of itinerant electrons.®®’®^ 

The above considerations suggest that the antisymmetric 
exchange couplings in Eq. (7) stabilize a specific spin con¬ 
figuration with a preference on the direction of magnetic mo¬ 
ments. Eor the ferromagnetic configuration in the same sub¬ 
lattice, the z-component order is preferred, which results in 
the up-down type antiferromagnetic state (z-UD) at and near 
half filling [see Eig. 2(a)], and the ferromagnetic state (z-E) at 
low and high fillings. On the other hand, for the antiferromag¬ 
netic configuration for the same sublattice, the jc-component 
order is favored, leading to the up-up-down-down type anti¬ 
ferromagnetic state (x-UUDD) irrespective of the sign in the 
effective exchange coupling between different sublattices [see 

Fig. 2(b)]. 

In the following sections, we focus on the z-UD order, 
since it accompanies an odd-parity multipole order composed 
of magnetic toroidal and quadrupole components.In this 
case, the ferroic toroidal order T is induced in the x direc¬ 
tion: T oc (V/Vpot) X (Si). In the following calculations, for 
simplicity, we restrict ourselves to the ordered states with the 
ordering vector, q — {q^, 0,0). 

3.2 Band deformation and magnetoelectric effect 

In this subsection, we examine the nature of the z-UD or¬ 
dered state in Eig. 2(a). As shown in the previous study, 
the z-UD order is expected to cause a band deformation with 
a shift of the band bottom. Also, the system may exhibit the 
magnetoelectric effect in which a staggered magnetic moment 
is induced by an electric current. We note that the magne¬ 
toelectric effect takes place even in the paramagnetic state. 
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due to the presence of the antisymmetric exchange couplings. 
We here confirm that these effects occur also in our extended 
Kondo lattice model in Eq. (7). 



k 


Fig. 3. Energy dispersion along the kx direction of the Hamiltonian in 
Eq. (7) under the z-UD order, as shown in the inset. Each band is doubly 
degenerate. The result is calculated at = \, = 0.1, ^3 = 0.2, = 0.2, 

7 = 4, and G = 0.5. The dotted red lines indicate the band bottoms. The 
dashed black curves represent the energy dispersion at G = 0. 


Figure 3 shows the band structure in our model in Eq. (7) 
along the kx direction from {-n, 0,0) to (tt, 0,0). The result is 
calculated by assuming the presence of the z-UD type antifer¬ 
romagnetic order with full polarization in localized electrons: 
Spi = (-l)'(z/2). We take fi = 1, f 2 = 0.1, fs = 0.2, U = 0.2, 
7 = 4, and G - 0.5. We also show the result at G = 0 for com¬ 
parison. As shown in Fig. 3, the z-UD order in the presence 
of the antisymmetric exchange couplings modifies the band 
structure in a peculiar way, with shifting the band bottom in 
the kx direction. This is similar to the previous results: 
although both spatial-inversion V and time-reversal T sym¬ 
metries are broken in this ordered state, the combined VT 
symmetry is retained, which ensures Ea-{k) — E-a-{k) (two¬ 
fold degeneracy in each band), while there is no guarantee to 
satisfy Ea-{k) - Ea-(-k). 

In our extended Kondo lattice model in Eq. (7), this band 
deformation is understood as follows. The Hamiltonian for 
the z-UD ordered state is represented by 

'Hchain=( M-v )’ ^^0) 

where u - S 2 + Gcrsf /4 and v = {Jcrll + V7Gsi:^)/2 
(cr = +1 for up and down spins); si = -2fi cos(kj/2) and 
S2 = -2f2 cos kx - 2f3 cos ky - 2f4 COS k^ represent the en¬ 
ergy dispersions for the different and same sublattices, respec¬ 
tively. We take the basis as (cAko-, cma), where CAka-ic^ka-) is 
the annihilation operator at A(B) sublattice with wave vector 
k and spin cr. By diagonalizing the Hamiltonian in Eq. (20), 


the eigenvalues are obtained as 

E(k) = £2 + ± (21) 

This indicates that the band deformation with the band bottom 
shift is due to VTGsj:^, which comes from the DiSk^ term in 
Eq. (7). Equation (21) also indicates that the peculiar band 
deformation does not occur for G = 0 (the standard Kondo 
lattice model), as shown in Fig. 3. 



He 

Fig. 4. Correlation between the z-UD moment and electric current in the x 
direction, K^, in Eq. (22). The data are obtained for several values of t[ and 
t 2 indicated in the figure, at ?3 = 0.2, (4 = 0.2, J = 4, G = 0.5, T = 0.01, 
and T] = 0.01. The inset shows a schematic picture for the magnetoelectric 
response expected from the correlation . 


Next, let US discuss the magnetoelectric effect by comput¬ 
ing the linear response of the staggered magnetization in the 
form of z-UD order by an electric current in the x direction. 
The response is calculated by the Kubo formula in the form: 

^(.s) _ J_ y /(g«t) - fiSmk) ^22) 

2 iVs ^ ^ ^nk ^mk ^nk ^mk W 

where is the system volume, /(e) is the Fermi dis¬ 
tribution function, cr"'” - {nA:|(-l)'cA|mA:), and 7"^ = 

{mk\Jx\nk} is the matrix element of the current operator Jx - 
(elh)&E(ex-KLMldkx', s„k and \nk) are the n-th eigenvalue and 
eigenstate of 77ex-KLM by assuming the z-UD order with full 
polarization in the localized spins. We set gji^ellh = 1. Fig- 

/o'l 

ure 4 shows the result of K„ as a function of the electron 
density tie - (l/N) c-^), where N is the total number 

of sites. The results are calculated at = 0.2, f 4 = 0.2, 7 = 4, 
G = 0.5, while changing f] and t 2 , which may correspond to a 
deformation of the zig-zag chain (see below). We set T = 0.01 
by taking the damping factor rj = 0.01. As shown in the inset 
of Fig. 4, the staggered moments in the z direction are induced 
by the electric current in the x direction, consistent with the 
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previous study.This is due to the presence of the toroidal 
component T ' in the multipole order. 

As shown in Fig. 4, the magnetoelectric response tends 
to be larger for larger fi/f 2 - This tendency is clearly seen in 
the low density region where the Fermi surface has a sim¬ 
ple shape. This suggests that the shallower zig-zag structure, 
where the intersublattice hopping becomes dominant, gives 
rise to a larger magnetoelectric effect. It also implies that 
the magnetoelectric response can be controlled by a uniax¬ 
ial pressure along or perpendicular to the zig-zag chain. The 
results are consistent with those in the single-band Hubbard 
model including the effect of the antisymmetric spin-orbit 
coupling. 

Moreover, the response also depends on J and G in a pe¬ 
culiar manner (not shown here). The G dependence is rather 
simple, at least, in the low density region. This is because the 
increase of G results in the increment of the antisymmetric 
spin-orbit coupling as in the third term of Eq. (7), which is the 
origin of the staggered magnetoelectric response. Meanwhile, 
the J dependence is complicated because the increase of J 
enhances not only the antisymmetric spin-orbit coupling but 
also the exchange coupling as in the second term of Eq. (7). 

3.3 Phase diagram 

In the previous subsection, we have studied the electronic 
and transport properties in the assumed z-UD ordered state. 
Now, we examine when and how such an ordered state is 
realized in the model in Eq. (7). Eor that purpose, we per¬ 
form three numerical calculations which are complementary 
to each other; variational calculation for the ground state in 
Sect. 3.3.1, simulated annealing in Sect. 3.3.2, and Monte 
Carlo simulation at finite temperatures in Sect. 3.3.3. 

3.3.1 Variational calculation 



2--4U4D(1) z-4U4D(2) 



J x-4U4D{l) x-4U4D(2) 


Fig. 5. Schematic pictures of ordering patterns of localized spins used in 
the variational calculations. U and D represent up- and down-spin polariza¬ 
tion, respectively, and F represents the ferromagnetic state. The prefices x and 
Z indicate the direction of the magnetic moments. 


First, we examine the ground state of the model in Eq. (7) 
by a variational calculation. Namely, for each parameter set 
of the Hamiltonian, we compare the zero-temperature grand 
potential per site, Q = E-^n^ {E - IN is. the inter¬ 

nal energy per site and fi is the chemical potential) for differ¬ 
ent magnetically-ordered states and determine the most stable 
one which gives the lowest Q. In the present calculations, we 
assume a collection of typical magnetic orders in the localized 
spins, up to the eight-site unit cell in the single chain as can¬ 
didates for the ground state; the states considered in this study 
are shown in Fig. 5. We consider only uniform ^ = 0 orders 
for all these patterns; namely, we assume that each magneti¬ 
cally ordered pattern is composed of a uniform arrangement 
of the magnetic unit cell in all the directions. The data are 
computed by approximating the integral over the folded Bril- 
louin zone using the sum over grid points of 64 x 64 x 64. 
The phase separated regions are determined by the method in 
Ref. 68. 

Figure 6 shows the phase diagram obtained by the varia¬ 
tional calculation, as a function of n^ and J. Figure 6(a) cor¬ 
responds to the result for the standard Kondo lattice model 
without the antisymmetric exchange couplings, i.e., at G = 0. 
In this case, there is no spin anisotropy because the model in 
Eq. (7) retains the spin rotational symmetry. In the low and 
high density regions, the ferromagnetic metallic phase ap¬ 
pears and becomes wider as J increases. The ferromagnetic 
phase is stabilized by the double-exchange mechanism.®’ On 
the other hand, a staggered antiferromagnetic order along the 
chain is stabilized at and near half-filling (We = 1) due to the 
effective antiferromagnetic interaction mentioned in Sect. 3.1. 
Note that an incommensurate order might take over the anti¬ 
ferromagnetic state in the weak-coupling region, which can¬ 
not be described in the present variational scheme within the 
limited sizes of magnetic unit cells; we will reexamine this 
point in the following subsections. In the intermediate n^ re¬ 
gion, there are several phases characterized by other order¬ 
ing wave vectors: UUDD and 4U4D antiferromagnetic orders 
(see Eig. 5). Note that, in all these phases, the common direc¬ 
tion of the magnetic moments can be taken arbitrarily owing 
to spin rotational symmetry of the system. 

Next, we discuss the effect of the antisymmetric exchange 
couplings by turning on G. When G ^ 0, the antisymmet¬ 
ric exchange couplings introduce the spin anisotropy, and 
hence, the system has a preference in the direction of the mag¬ 
netic moments in each phase. Figure 6(b) shows the result at 
G = 0.5. At and near half filling, the quantized axis in the 
up-down antiferromagnetic phase prefers to be fixed in the z 
direction. This is the z-UD phase with multipole ordering dis¬ 
cussed in Sect. 3.2, which gives rise to a band deformation 
with a band bottom shift and exhibits the magnetoelectric ef¬ 
fect. The result in Fig. 6 indicates that this UD phase tends to 
be stabilized by G, as was discussed in Sect. 3.1. 

Meanwhile, the UUDD phases near quarter and three quar¬ 
ter fillings are also present at G = 0.5, as shown in Fig. 6(b). In 
this case, however, the magnetic moments are aligned within 
the xy plane by the spin anisotropy; we denote it by UUDD 
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Fig. 6. Ground-state phase diagram of the model in Eq. (7) obtained by the 
variational calculation at (a) G = 0 and (b) G = 0.5. Other pai'ameters ai‘e 
taken to be t\ = \, t 2 = 0.1, ^3 = 0.2, and ^4 = 0.2. PS indicates a phase 
separated region. Ordering patterns are shown in Fig. 5. 


with the prehx x, while the moment can point to any direction 
within the xy plane. On the other hand, in the ferromagnetic 
phases in the low and high hlling regions, the moments are 
polarized in the z direction. 

For all the phases we obtained, the direction of the mag¬ 
netic moments as well as the stable parameter region is con¬ 
sistent with the arguments for the two-site problem discussed 
in Sect. 3.1. The tendency does not change while changing 
tift 2 (not shown). 


3.3.2 Simulated annealing 

In order to conhrm the stability of the z-UD state by a 
more unbiased method than the variational calculation in 
Sect. 3.3.1, we adopt simulated annealing. The simulated an¬ 
nealing is an optimization method for hnding the global min¬ 
imum of a function that possesses many local minima.By 
using the method, we obtain the accurate magnetic order in 
the ground state within the unit cell we set in the calcula¬ 
tion; we do not need to assume a specihc magnetic order, in 
contrast to the variational calculation in Sect. 3.3.1. We opti¬ 
mize the classical localized spins by the simulated annealing 
as follows; temperature T is decreased gradually, and for each 
T, the spin conhguration is updated by the Monte Carlo sam¬ 
pling with the single-spin flip algorithm. 

In the present calculations, we performed the simulated an¬ 
nealing while decreasing T in a geometrical way, T„+i = aTn, 
where Tn is the temperature in the n-th step. We started from 
the initial temperature Tq - 1.0 by taking the coefficient of 
geometrical cooling a - 0.97 and the total steps of cooling 
270; the final temperature Trio reaches down to ~ 3.0 x 10“"^. 
We considered the systems with N - 16x1x1 and 24 x 1 x 1 
while introducing a supercell consisting of Nt - 8^ copies of 
the A^-site lattice to reduce the finite-size effect. 



Q.X 



Fig. 7. qx dependence of the spin structure factor at = 0 divided by 

the system size N obtained by the simulated annealing. See Ref.70. The data 
are calculated at = 1, ?2 = 0.1, ^3 = 0.2, 4 = 0.2, and 7 = 4. The results 
for G = 0 and G = 0.5 with the system size N = \6 and 24 are shown in each 
figure: (a) /i = -0.4 and (b) fj. = -0.8. Electron fillings for each /i are also 
shown in the figure. 
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Figure 7 shows the spin structure factor obtained by the 
simulated annealing. The spin structure factor is defined by 

ij 

where ^ is a wave vector and r, is the position vector at site i. 
In Fig. 7, S q/N is plotted as a function of - (qx, 0,0). 

Note that, for the perfect UD order, SqjN - 0.25 at g = 
{n, 0,0), and otherwise zero.™ At/j.- 0.0 (almost half filling, 
He - 1.0), the spin structure factor shows a sharp peak at g = 
(tt, 0,0) for both G - 0 and G - 0.5 (not shown here). The 
peak values do not depend on the system size, which suggests 
the UD order is stable even when allowing other magnetic 
orders in larger unit cells than in the variational calculations. 
Similar to the results in Sect. 3.3.1, the moments for the UD 
order are aligned in the z direction when introducing G. 

As shown in Fig. 7(a), the peaks remain at g = (tt, 0,0) 
against a slight decrease of Wg by a few percent. However, the 
peak for G - 0 slightly decreases (~ 1.5%) when increasing 
the system size from N = 16 to 24, while the peak value for 
G = 0.5 remains almost unchanged (the change is less than 
0.3%). These results imply that the UD state survives against 
a slight doping, and the antisymmetric exchange couplings 
stabilize the z-UD state. 

While further decreasing Wg, as shown in Fig. 7(b), the peak 
position for G = 0 shifts to a smaller qx. This suggests the 
commensurate UD state is no longer stable and taken over by 
an incommensurate order with a longer period, as anticipated 
in the variational arguments in the previous subsection. Mean¬ 
while, for G = 0.5, the z-UD state remains stable, as shown in 
Fig. 7(b). Thus, the results of the simulated annealing confirm 
the variational results in the previous subsection: the z-UD or¬ 
dered state appears as a stable phase at and near half filling in 
the presence of the antisymmetric exchange couplings. 

3.3.3 Monte Carlo simulation 

In this subsection, we investigate the stability of the z-UD 
phase at finite temperatures by Monte Carlo simulation. In the 
Monte Carlo calculations, we adopt the standard method for 
the spin-charge coupled systems with classical localized mo¬ 
ments.^* We typically performed 10,000-90,000 Monte Carlo 
steps after 10,000 steps for thermalization. The statistical er¬ 
rors were estimated by dividing the data into five to ten bins 
and calculating the standard deviation among the bins. The 
calculations were performed on the N - 4LxLx L-site lattice 
with L = 2 and 4 under the periodic boundary conditions in all 
the directions. For L = 2, we introduced a supercell consisting 
of Nk - 2? copies of the A-site lattice so that the finite-size 
effect is reduced. In order to stabilize the qy - q^ - 0 or¬ 
der as in the previous sections, we introduce the additional 
ferromagnetic exchange interaction between localized spins 
’JLp = ■ ^j''’ we take Jp - 0.1 and the sum of 

(h j)yz over the nearest-neighbor sites for the y and z direc¬ 
tions. 

Figure 8(a) shows the temperature dependence of the or- 



0.02 0.06 0.10 
Temperature 

Fig. 8. Monte Carlo results for (a) the order parameter for the ^-UD order, 
mg, and (b) the electron density at = -0.1 and fj. = -0.5. The calcula¬ 
tions were done at ?i = 1, ^2 = 0-C ?3 = 0.2, ?4 = 0.2 7 = 4, and G = 0.5 for 
the system sizes = 8 x 2 x 2 and 16x4x4. 


der parameter for the UD order, mg. Here, mg is defined by 
mg = [5g/A]*™, where 5g is the spin structure factor in 
Eq. (23) 24 Q - (tt, 0,0). The data are calculated for the two 
different system sizes at 7 = 4 and G = 0.5 for two different 
values of the chemical potential, p = -0.1 and p = -0.5. The 
results show that mg develops rapidly below a particular tem¬ 
perature. We confirmed that the magnetic moments are along 
the z direction by analyzing the spin component of S g. These 
indicate a phase transition from the high-temperature param¬ 
agnetic state to the low-temperature z-UD ordered state. The 
rough estimates of the critical temperature Tc can be obtained 
from the inflection point of mQ(T): Tc - 0.09 2 A q - -0.1 
(we ~ 0.96) and Tc - 0.075 24 q - -0.5 (n^ ~ 0.85) [see 
also Fig. 8(b)]. With further decreasing temperature, the or¬ 
der parameter mg approaches its saturated value 0.5 in the 
ground state. This confirms that the multipole ordered state 
found in the variational calculations remains stable against 
thermal fluctuations as well as the carrier doping. 

By calculating mg(r) while varying p in a similar way, we 
obtain the finite-temperature phase diagram in Fig. 9. The z- 
UD phase is stable around half filling, with a dome-like shape 
of Tc with its maximum around n^ - 1. Although it is diffi¬ 
cult to determine the phase boundary at n^ ~ 0.7 within the 
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Fig. 9. Finite-temperature phase diagram of the model in Eq. (7) as a func¬ 
tion of the electron density Hq obtained by the Monte Carlo simulation. The 
symbols indicate the critical temperatures, which are estimated by the inflec¬ 
tion points of the order parameter mg plotted in Fig. 8. The horizontal eiTor 
bars indicate the statistical en'ors of rie at the chemical potential for the crit¬ 
ical temperature. The results are calculated at ?! = I, t 2 = 0.1, ts = 0.2, 
t 4 = 0.2, 7 = 4, and G = 0.5. 


limited system sizes in the current calculations, the result indi¬ 
cates that the multipole ordered state is robustly stable around 
half filling, consistent with the results by the variational cal¬ 
culation and simulated annealing. 

4. Summary 

In summary, we have investigated the odd-parity multipole 
ordering that is spontaneously induced by the antisymmetric 
spin-orbit coupling in the systems with local parity mixing. 
Starting from a general form of the site-dependent antisym¬ 
metric hybridization between different parity orbitals, we de¬ 
rived an effective low-energy model with site-dependent anti¬ 
symmetric exchange couplings. This is an extended Kondo 
lattice model, which is a fundamental model for consider¬ 
ing the effect of antisymmetric hybridization in d- and /- 
electron systems. We have analyzed the model on a quasi-one- 
dimensional zig-zag lattice, as the minimal lattice structure 
describing local parity mixing. We found that the antisymmet¬ 
ric exchange couplings induce the effective hopping of con¬ 
duction electrons depending on the configurations of localized 
spins. One of the stable configurations is a Neel type antifer¬ 
romagnetic order with the moments perpendicular to the zig¬ 
zag plane. The Neel order accompanies an odd-parity mul¬ 
tipole order composed of magnetic toroidal and quadrupole 
components. This unusual multipole order exhibits a band de¬ 
formation with a band bottom shift and magnetoelectric re¬ 
sponse, due to the activated toroidal moment. We have inves¬ 
tigated the stability of the multipole ordered state by the com¬ 
plementary numerical calculations, i.e., the variational calcu¬ 
lation for the ground state, the simulated annealing, and the 


Monte Carlo simulation at finite temperatures. We found that 
the multipole ordered state is indeed stabilized by the anti¬ 
symmetric exchange couplings in a wide parameter range at 
and near half filing. 

Our results will stimulate further studies of odd-parity mul¬ 
tipole ordering in systems with local parity mixing. Multipole 
orders similar to that in the present study will be widely ob¬ 
served in the materials to meet the following conditions: (i) 
local inversion symmetry breaking due to the lattice struc¬ 
ture (zig-zag, honeycomb, diamond,...), (ii) strong spin-orbit 
coupling, (iii) hybridization between different parity orbitals, 
e.g., s-f, p-d, and d-f. There are many candidate materials 
to meet such conditions in /-electron systems, as mentioned 
in the introduction. We anticipate a similar situation also in 
Ad- and StZ-electron systems where localized levels are ex¬ 
pected under some crystal field splitting. It is desired to sys¬ 
tematically study such systems from the viewpoint of odd- 
parity multipole ordering for further understanding of the ex¬ 
otic magnetism, accompanying peculiar electronic and trans¬ 
port properties. 

SH is supported by Grant-in-Aid for ISPS Fellows. This 
work was supported by Grants-in-Aid for Scientific Research 
(No. 24340076), the Strategic Programs for Innovative Re¬ 
search (SPIRE), MEXT, and the Computational Materials 
Science Initiative (CMSI), Japan. 

Appendix: Canonical Transformation in the Presence of 
Antisymmetric Hybridization 

In this Appendix, we derive the extended Kondo lattice 
model in Eq. (7) from the extended periodic Anderson model 
with the antisymmetric hybridization in Eq. (2). Eollowing the 
procedure for deriving the standard Kondo lattice model by 
the Schrieffer-Wolff transformation,®^ we treat the hybridiza¬ 
tion Tfi in Eq. (4) as a perturbation to 7/) in Eq. (3), and 
perform the second-order perturbation expansion by using 
a canonical transformation. The canonical transformation is 
represented by 

(Al) 

where S is the anti-Hermitian operator, determined so as to 
satisfy the following relation: 

7/i-t [ 5 , 7 / 0 ] = 0. (A-2) 

Here, [■ ■ ■ ] represents a commutator. Expanding Eq. (AT) up 
to the second-order in 7/i, we end up with an effective model 
of Kondo lattice type: 

'7/ex-KLM=7/o + i[5,7/l]. (A-3) 

Note that S is Oi'Hi) due to the relation in Eq. (A-2). 

Specifically, the operator S is obtained in the form: 

5= ^ [Ai^Ak)Bpv.Ak)cljpv^,-'^-c], (A-4) 
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where 


Alik) = 


Aiii(k) Al^(k) j ’ 


= Viik)cr^ + gf(k) ■ 0-, 
BpUk) = {p(*) + 


Pik) = 
Qik)^ 


1 


eik) - £o ’ 
1 


1 


eik) - Eq eik) - Eq-U 


(A-5) 

(A-6) 

(A-7) 

(A-8) 

(A-9) 


fi(^) = 2//' ^ii'ik) is the energy dispersion of conduction elec¬ 
trons, which is obtained by the Fourier transform of the first 
term in Eq. (2) [see below in Eq. (A-10)], and Rp is the posi¬ 
tion vector for unit cell p. By substituting S in Eq. (A4) into 
the second term in Eq. (A-3), we obtain the effective Hamil¬ 
tonian, which is represented by 


"ktex-KLU - 


Z ‘^U'ik)cl^Cp,^ + ^ Jlk'kSpl ■ Siic'k^R^ 

l,p,k,k' 


terms for / electrons proportional to 

Equation (A-10) provides a general form of the extended 
Kondo lattice model applicable to any lattice structure. Ei- 
nally, we simplify the k dependence in Eq. (A-10), bearing 
a quasi-one-dimensional system composed of zig-zag chains 
[see Eig. 1(b)] in mind. Namely, we drop all the k depen¬ 
dences except for the most interesting one, sin kx in ik) in 
Eq. (6), as 

Jlk'k —> (A-14) 

^ (A-15) 

^m'k ^ (A-16) 

where 6pv is the Kronecker delta. Here, J and G are gener¬ 
ally positive from Eqs. (A-11) and (A-13). Note that Di de¬ 
pends on sublattice index, which comes from g'^/ik) oc (-1)^ 
for the present zig-zag lattice structure. We also approximate 
the magnitude of Di by V7G from Eqs. (A-11), (A-12), and 
(A-13). By substituting Eqs. (A-14), (A-15), and (A-16) into 
Eq. (A-10), we obtain the extended Kondo lattice Hamilto¬ 
nian for the zig-zag lattice structure as in Eq. (7). 


X [{^^k'k pl’^lk'ki + Spitlik'k^ + 2i-5p/-4't) 


l,p,k,k' 


l,p,k,k’ 


l,p,k,k' 



[S],inik'ki -Spitiik'k^ + 2SZsZ'k) 

1) 

2) 


^Z^k’k ~ ^pl^ik’k + ^pl^lk'k)]SRp + H.C.j 

3) 



4) 

[Gi^S 

pl4k'k + ^pl^lk'k ~ 

5) 


^Z^^k'k + ^ pl4k'k + 

6) 


S;iSik'k + S-p,sZ,,-2SZsZ,,)}6R, 

7) 

[{-GTk'k 

[S^nik’ki - Spinik’kt - 2S^^isZk) 

8) 


{^pl^ik’k ~ ^pl^Ik'k + ^Z^lk'k) 

9) 


{shmk'ki + Spiiiik'k^ - 2iSl,sZ,^)]6R^ + H.C.] 

10) 

’ll) 


(A-IO) 


where — 2 a-tvith “ ^Ik'o'^lkcr^ ^Ik'k — ^^Ik'k"^ 
sJk,k)l'2, and = (^tk’k ~ %'t^/2i- The coefficients for the 
exchange couplings are given by 

Jlk'k - yiik')V*iik)Qk,k', (A-11) 

Dik'k = Viik')[gfik)rQk.k', (A-12) 

= [gf^k')][gf\k)rQk,k', (A-13) 


^Ik'k 


where Qk,k' - Qik)+Qik’) and p,v = x, y, z. In the derivation, 
we used the condition = 1, assuming the / electrons are 
well localized in a singly-occupied state at each site, as the 
standard Kondo lattice model. In the same spirit, we neglect 
pair hopping terms proportional to and one-body 
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